Mechanism for puddle formation in graphene 
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When graphene is close to charge neutrality, its energy landscape is highly inhomogeneous, forming 
a sea of electron-like and hole-like puddles, which determine the properties of graphene at low carrier 
density. However, the details of the puddle formation have remained elusive. We demonstrate 
numerically that in sharp contrast to monolayer graphene, the normalized autocorrelation function 
for the puddle landscape in bilayer graphene depends only on the distance between the graphene 
and the source of the long-ranged impurity potential. By comparing with available experimental 
data, we find quantitative evidence for the implied differences in scanning tunneling microscopy 
measurements of electron and hole puddles for monolayer and bilayer graphene in nominally the 
same disorder potential. 
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I. INTRODUCTION 

In monolayer graphene, the hexagonal arrangement of 
carbon atoms dictates that in the absence of atomic-scale 
disorder, graphene is a gapless semiconductor—!^ that is 
always metallic at low temperature,- This metallic be- 
havior holds even in the presence of quantum interfer- 
ence and strong disorder,** in stark contrast to most other 
materials, which undergo a metal-to-insulator transition 
at low carrier densityi^i^ The physical origin for this ro- 
bust metallic state is that the ground-state of graphene 
at vanishing mean carrier density becomes spatially in- 
homogeneous, breaking up into electron-rich and hole- 
rich metallic regions connected by highly conducting p- 
n j unctions Bilayer graphene comprising two sheets of 
graphene that become strongly coupled due to the AB 
stacking arrangement^ shares some properties with reg- 
ular semiconductors (such as the parabolic band disper- 
sion) and in other ways behaves like monolayer graphene, 
including having chiral wavefunctions and forming elec- 
tron and hole puddles at low density. 

These electron and hole puddles have now been ob- 
served in several experiments of exfoliated graphene on 
an insulating Si02 substrate including Refs.[9l-[l4l While 
these authors suggest that long-range charged impurities 
in the substrate could be responsible for the spatial inho- 
mogeneity, detailed comparisons to microscopic models 
have not been made. 

In this paper, we demonstrate that differences between 
the spatial properties of puddles in monolayer and bi- 
layer graphene can be quantitatively explained by the 
differences in the screening properties of the two systems 
(that ultimately arises from the differences in their band- 
structure). Numerical results show that the correlation 
length for bilayer graphene is relatively independent of 
density and significantly smaller than that of monolayer 
graphene for a typical range of impurity densities. Fi- 
nally, we find good quantitative agreement when com- 
paring our results with available experimental data. 



The rest of the paper is organized as follows. In SecHIl 
we outline the the theoretical model, providing a heuris- 
tic understanding of our results using the Thomas-Fermi 
(TF) screening theory. However, the TF significantly 
underestimates the effect of electronic screening in both 
monolayer and bilayer graphene. It is therefore neces- 
sary to use the Random Phase Approximation (RPA) 
screening theory, which we discuss in Sec. IIIIl Our main 
finding is that the puddle correlation length in bilayer 
graphene (^ w 3.5 nm) is relatively insensitive to the 
impurity concentration and carrier doping. This is in 
contrast to monolayer graphene, where the puddle corre- 
lation length varies from 3 nm in dirty samples to more 
than 35 nm in clean samples. 

The comparison with experiment is done in Sec. IIVI 
where we examine three different experimental results: 
(1) We consider first the experimentally determined nor- 
malized correlation function A{r) (see definition below) 
obtained from the scanning tunneling microscopy (STM) 
data reported for exfoliated bilayer graphene in Ref. [ij. 
The full functional form of A{r) agrees with the theory 
where the only adjusted parameter in the theory is the 
distance d of the impurities from the graphene sheet. 
In particular the experimentally determined correlation 
length ^ = (3.68±0.03) nm, defined here as the half-width 
at half-maximum (HWHM) decay length of A{r) agrees 
well with the d = 1 nm RPA theory value of ^ = 3.4 nm. 
This value of d is both reasonable and consistent with 
those determined from other transport measurements on 
bilayer graphencjii^ (2) From the mon olay er graphene 
STM experimental data reported in Ref. [Tj, we extract a 
correlation length 6 nm < ^ < 11 nm. Since the measure- 
ment a^'^d^ were made on the same exfoliated graphene 
sample containing both single layer and bilayer graphene 
regions, we expect that the extrinsic disorder potential 
is statistically identical for the two samples. Therefore, 
using the value of d = 1 nm (discussed above) and the 
disorder induced Dirac point shift reported in Ref. [l3, 
we calculate theoretically (without any adjustable pa- 
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rameters) that a monolayer graphene sample in the same 
disorder environment would have a puddle correlation 
length ^ = 8 nm, in reasonable agreement with the exper- 
iment. (3) We then compare A{r) obtained using scan- 
ning Coulomb blockade spectroscopy reported in Ref. [12 
with our theoretical results for monolayer graphene at 
the Dirac point. The parameters used in the theory 
were obtained from separate transport measurements on 
the same experimental sampleJ^ The agreement between 
theory and experiment is remarkable since it involves no 
adjustable parameters. Finally in Sec. |Vl we conclude 
by making predictions for future experiments involving 
monolayer and bilayer graphene on BN substrates. 



II. FORMALISM 

The doping level of graphene can be measured in a va- 
riety of different ways. In transport measurements, the 
gate voltage potential that yields the resistivity maxi- 
mum identifies the extrinsic doping level due to extrane- 
ous sources, such as charged impurities in the substrate 
impurities with density, Uimp- While the width of the 
resistivity maximum is a measure of the homogeneity of 
the sample,-- or the electron-hole puddle distribution. 

In local probe measurements, such as STM, the Dirac 
point energy relative to the Fermi-level can be observed 
as a minimum in the tunneling differential conductance, 
dl/dV as a function of tunneling bias. Knowing the 
electronic dispersion relation (see details below), for a 
particular gate voltage Vg, this spatial map V{r) of the 
Dirac point variation can then be used to extract the spa- 
tial distribution of the local carrier density (characterized 
by a width nmis)- We can also characterize the puddles 
through the radially averaged autocorrelation function 

27r 
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where the angular brackets denote an average over the 
image area and the ^-integration averages over orienta- 
tions. 

Notice that while C(0) = T/^^g (which is related to 
rirms) characterizes the fluctuations in the puddle depth, 
C{r) describes the spatial profile of the electron and hole 
puddles. We find it useful to consider the normalized cor- 
relation function A{r) — C{r)/C{0). We will argue below 
that A{r) and C(0) are quite different physical quantities 
that depends quite differently on the parameters of the 
extrinsic impurity potential. (In addition, for a typical 
STM experiment, where the shifts in the Dirac point are 
determined^*^ from the shifts in dl/dV at fixed Vg, the 
determination of C (0) is complicated by the experimen- 
tal uncertainty in converting spatial maps of dl/dV to 
Dirac point energy shifts V{r). By contrast, the much 
smaller uncertainty in A{r) is mostly determined by the 
spatial resolution and image area.) 



To theoretically compute the correlation functions for 
puddles in graphene, we make two assumptions. First, 
the impurity potential comes from a random two dimen- 
sional distribution of charged impurities displaced by a 
distance d from the plane with density nimp. This model 
has been highly successful in describing the effect of disor- 
der in semiconductor heterojunctionsi^ and in graphencii 
This spatially varying potential gives rise to a varying 
charge density and local variations in the screening of 
the potential. Second, we assume that it is possible to 
find a global screening function e(g, rieff) that adequately 
describes the effects of these local screening variations. 
Here, the screening depends on the disorder potential 
only through an effective carrier density rics- This self- 
consistent screening model has been used previously to 
understand the minimum conductivity problem in both 
monolayer— and bilayer graphene.— If these two assump- 
tions are satisfied, the correlation function, aside from a 
prcfactor of rtimp, then depends only on the screened im- 
purity potential 

r( ^ 9 ^'V gexp(-2gd) 

C(r) = 27rnimp — / dq— Jo(gr),(2) 

where k is the bulk (3D) dielectric constant, e{q) is the 
surface (2D) screening function in the plane, — e the elec- 
tron charge, and Jo{x) is a Bessel function. 

As an illustration, consider the Thomas-Fermi (TF) 
screening for which the surface dielectric function is given 
by e(g,ncff) = 1 + (?TF("-off)/g, where gTF("-off) (dis- 
cussed below) is the Thomas-Fermi screening wavevec- 
tor. Shown in Fig. [T] is a calculation of C(r) for different 
values of gxF and d. Notice that for fixed ^tf, the spa- 
tial dependence of C(r) depends on d and qtf, but not 
on riimp- On the other-hand, the function C(0) (which, 
within the TF can be calculated analytically) depends on 
JT-imp, K, d and gxF- This is why we find it useful to use 
the normalized correlation function A{r) — C(r)/C(0) 
that describes the spatial profile of the screened impu- 
rity potential. We also emphasize that A{r) contains 
different information than the typical puddle size - for 
example, the puddle correlation length ^ (recall that ^ is 
defined as the HWHM of A{r)) describes the width of the 
screened impurity potential, and not the mean impurity 
separation. For example, in the low impurity density 
limit, spatial maps of the puddles would show isolated 
impurities, but A(r) would not change (for fixed (?tf)- 
Moreover, since A{r) scales differently with c^tf and d, 
in principle, both of these length scales can be extracted 
from a measurement of A{r). 

Within the TF screening theory, any differences be- 
tween monolayer and bilayer graphene can only arise 
from differences in (jTF("-cff)- The linear dispersion in 
monolayer graphene and the hyperbolic dispersion in bi- 
layer graphene gives rise to these differences. The low 
energy linearly dispersing bands of monolayer graphene 
can be modeled by a single parameter, the Fermi ve- 
locity Up, or equivalently, the effective fine-structure 
constant = e'^/{Khv-p) sa 0.8. characterizes the 
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FIG. 1: Theoretical calculations for the correlation function 
C(r) using the Thomas- Fermi approximation. This autocor- 
relation function depends separately on the typical distance 
d of long-ranged impurities from the graphene sheet, and ^tf 
the inverse effective screening length, allowing them to be de- 
termined independently. Symbols show the correlation length 
^, defined as the IfWHM length. 
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FIG. 2: Effective carrier density as a function of impurity den- 
sity assuming d — I nm, rs = 0.8, and no = 2.3 x 10^^ cm~^. 
For bilayer graphene, the blue circles show the Thomas-Fermi 
approximation and the red squares are RPA results. The em- 
pirical relation tIch = -yrtimp ni adequately captures the RPA 
results, with ni = 6.8 x 10^^ cm"'^. 



strength of the electron-electron interaction for graphene 
on a Si02 substratei and is useful because we are in- 
terested in the screening properties of graphene. The 
Thomas-Fermi screening wavevector is related to the den- 
sity of states, and for monolayer graphene is given by 
iTpincg) = Argy/Trneg, where is the effective carrier 
density. 

Bilayer graphene can be modeled with a hyperbolic 
dispersion with two parameters v-p (throughout this 
manuscript, vp is the Fermi velocity of a single decoupled 
graphene sheet), and the low-energy effective mass moff. 
For simplicity we use for the two parameters (defined 
above) and uq = mlgV^/{h'^TT) w 2.3 x 10"^^ cm~^ which 
is the characteristic density scale for the crossover from a 
(low density) parabolic to a (high density) linear disper- 
sion. The Thomas- Fermi screening wavevect or for bilaye r 
graphene is given by gTF(n-cff) = ^rsy/irnoy^l + n^s/nQ. 

As the system approaches the Dirac point, the fluc- 
tuations in carrier density become larger than the av- 
erage density. In this case, screening varies spatially 
with the density fluctuations. We assume that it is pos- 
sible to describe the effect of this screening by using 
the screening for an ideal system and an effective car- 
rier density ricff obtained self-consistentlyJ^ This is done 
by equating the squared Fermi level shift with respect 
to the Dirac point with the square of the potential fluc- 
tuations, E'^[n = Ticff] = C'(O) where C(0) is defined in 
Eq. [2]and E[ri\ — hvp^/wn for monolayer graphene, and 

for bilayer graphene. 



E[n] = v|moff a/1 -I- ti/tiq - 1 
The result of this procedure are shown in Fig. [21 

Within the TF theory, we can now qualitatively dis- 
cuss the main differences between monolayer and bilayer 
graphene. For bilayer graphene, the inverse screening 
length changes only slightly from the low density value of 
9BLG ~ 4rs->/7rno that is set entirely by the band param- 



eters. As a consequence, the puddle correlation length 
does not change with the impurity concentration or car- 
rier density, and depends only on the distance d of the 
bilayer graphene sheet from the source of the long-ranged 
impurity potential. In contrast, for monolayer graphene, 
giviLG = '^TsyJirricG depends essentially on C(0) (and 
therefore on nimp). This heuristic description (which we 
make more quantitative below) implies that depending 
on the sample quality, choice of substrate, or doping, the 
puddle correlation length in monolayer graphene (but not 
bilayer graphene) could vary by more than an order of 
magnitude. 



III. RANDOM PHASE APPROXIMATION 

While the TF screening theory discussed in the previ- 
ous section is useful to obtain a qualitative picture, we 
find that it significantly underestimates the effect of elec- 
tronic screening. In both monolayer and bilayer graphene 
it gives larger values for C(0) and smaller values for ^. 
In what follows we use the the random phase approxi- 
mation (RPA) where the screening function is obtained 
using e(g) = H- qTF\l[q)/q. 

The normalized polarizability n(g) for monolayer— 
and bilayeri'^S graphene are both available in the litera- 
ture. We note that for monolayer graphene n((j) depends 
only on the dimensionless variable x — q/(2^7rncff ) 



n(a;) = l + e'(a;-l) 



TTX X , , \J\ — X 

arccsc X 

4 2 ^ ' 2 



B(l-x)^B{x-\) — , 



(3) 



where Q(x) is a step- function and -Do is the density of 



4 



states with qtf = 2Tr{e'^ /k)Do = i^irncers- 

In contrast, for bilayer graphene, the polarizabihty de- 
pends both on the scaled momentum transfer x, and on 
V = iT-cff/nQ < 8, which parameterizes the bilayer hyper- 
bolic dispersion relation. For 77 ^ 1, the bilayer graphene 
dispersion is quadratic, while for 1 ^ ?7 < 8, the disper- 
sion is linear. For > 8, one must consider the effects 
of a second higher-energy band that provides additional 



screening-^ and is not considered here. By restricting 
the density to n < 8noj can simplify the expression 
for the bilayer polarizabihty reported in Ref. Us- 
ing the bilayer density of states D — Dq^/I + rj, with 
Dq — 2m/(7r/i^), the normalized polarizabihty function 
n(a;,77) = Il{q)/DQ is given by 



n(x,77) 



f{x,ri) + 0{x - l)g{x,ri), 



l^x^r] + x'^{2 + r] + 2^1 + ri) 



1/2 



In 
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r]x^ 



-1 + VT+Tj 



1 3X^T] - 1 

- H 

2 2x^ 



arctan(a;Y^) 



+ \/l — x'^rj ^2arctanh(v/l — ryx^) — arcsinh 
-Va;2 - 1(1 + 7] + 2x^7] - y/TTlj) ^ {3x^v - 1) 

2Xy^ 
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(4) 



The polarizabihty functions for monolayer and bilayer 
graphene are shown in Fig. [31 What is left is to cal- 
culate the effective residual density rieff as a function 
of impurity concentration. As discussed earlier, this is 
obtained by first calculating the autocorrelation func- 
tion C(0) from Eq. [2J using the RPA results shown in 
Fig. [3] Figure |4] shows the autocorrelation function C(0) 
obtained using the RPA results (solid lines) as well the 
Thomas- Fermi results (dashed lines). We note that ex- 
cept at very high density (where both monolayer and bi- 
layer graphene approach the "complete screening" limit, 
with C(0) = (4fcFrsd)~2), the Thomas-Fermi approxima- 
tion grossly underestimates the effect of screening. More- 
over, for typical densities in bilayer graphene C(0) is ap- 
proximately constant and independent of carrier density 
consistent with the heuristic picture discussed at the end 

of Sec. nn 

The effective density calculated within the RPA is 
shown in Fig. [5J For bilayer graphene, we find that the 
following empirical relationship adequately describes the 
numerical results 

n^s = ^ninip ni, (5) 

where n\ = 11.5 x 10^^ cm"^ for the Thomas-Fermi ap- 
proximation, and Til = 6.8 x 10^^ cm^^ for the RPA 
results. The scaling of the bilayer effective density rioff ^ 
^riiijip can be anticipated for the TF approximation in 
the limit n^s ^ no. However, it is surprising that the 
simple empirical relation continues to hold both for the 
RPA screening theory, and for larger values of noff . 
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FIG. 3: The Random Phase Approximation polarizability 
function 11(g) normalized by the density of states for mono- 
layer graphene and for bilayer graphene with a hyperbolic dis- 
persion. Also shown is the parabolic approximation for the 
bilayer, which can be obtained from the hyperbolic dispersion 
when T] <^\. The Thomas- Fermi approximation discussed in 
the text corresponds to the assumption that the normalized 
n(g) = 1 for all g. 



This dependence of ricg ~ ^nimp for bilayer graphene 
should be contrasted with similar results obtained 
previously for monolayer graphene where rioff — 
2r^C[0]7T,imp cannot be captured by a similar empirical 
fit. For comparison, these earlier results are also shown in 
Fig. [21 where we emphasize that for a given impurity con- 
centration (riimp), bilayer graphene exhibits larger den- 
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FIG. 4: Potential autocorrelation function C[0] for mono- 
layer and bilayer graphene. At the Dirac point, ftp is the 
Fermi wavevector arising from the effective carrier density 
i.e. fcp = y/nrics. For large density Akprsd ^ 1, both the 
monolayer and bilayer results approach the "complete screen- 
ing" limit, defined here as C(0) = (4fcFrsd)~^. Notice that 
the Thomas-Fermi approximation shown as dashed lines cap- 
tures the correct qualitative behavior, but can give signifi- 
cantly larger values for C[0], and is therefore unsuitable for 
quantitative comparisons. 



sity fluctuations (fics) than monolayer graphene. Finally, 
using Eq. [21 we can also calculate the puddle correla- 
tion function, and the corresponding HWHM correlation 
length, ^, that is shown in Fig. [5] 
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FIG. 5: Theoretical results for the puddle correlation length 
at the Dirac point as a function of impurity concentration. 
While the puddle size in bilayer graphene « 3.5 nm) is 
relatively insensitive to the disorder concentration, the size of 
the puddles in monolayer graphene varies from 3 nm in dirty 
samples to more than 35 nm in clean samples. 



IV. COMPARISON WITH EXPERIMENTS 

We now compare the calculated correlation functions 
with experiment. Figure |6l^a) shows that for bilayer 
graphene, A{r) extracted from the data reported in 
Ref. [13 agrees with the calculation for d = 1 nm.^^- The 
circles show the experimental data and the RPA theory 
for bilayer graphene is shown for d = 1 nm (solid curve) 
and d — 0.5 nm (dashed curve). The theoretical results 
are insensitive to the impurity concentration Tiimp and to 
how far the doping is away from the Dirac point. Con- 
sequently, the only free parameter in the theory is the 
distance d of the impurities from the graphene sheet. 

In Ref. [3, we reported a maximum peak-to-peak car- 
rier density fluctuation of 3.6 x 10^^ cm~^. To extract an 
impurity density (using the results in Fig. we need to 
estimate rii-ms from this peak-to-peak value. By assuming 
that the carrier density has a Gaussian distribution, and 
estimating that the peak-to-peak corresponds to a mea- 
surement of 4cr, we roughly estimate that rics — Jirms ~ 
10^^ cm~^ and that riimp = n^^^/ni « 1.2 x 10^° cm^^ 
for the substrate induced impurities in that experiment. 

We can use this value to predict theoretically the corre- 
sponding density fluctuations in the adjacent monolayer 
sample reported in Ref. 13. However, one complication 
is that the theory discussed in Sec. |in]was developed for 
monolayer graphene at the Dirac point, while the exper- 
imental data was taken at a backgate induced density 
Ug — 1.4 X 10^^ cm~^. Very far from the Dirac point, 
i.e. when riimp/ng — ?> 0, the potential fluctuations Vrms 
can be obtained from Eq. [2] by setting rics = rig on the 
right-hand side. In this case, the density fluctuations are 



(6) 



We note that when , 



^ 1, we can use the result 
(see Fig. 131) to obtain rii-^s ~ \/'T-imp/(87rc?2). 
However, these constraints are not fully satisfled in the 
experimental data. Calculating Tirms in the crossover be- 
tween the limits Ug = and Ug ^ nimp is more compli- 
cated. For our purposes, it is sufficient to extrapolate be- 
tween the low-density and high-density limits by adding 
the two contributions in quadrature, and solving for n^rns 
self-consistently. This procedure gives 



= 2r,Jni„,pCRPA(o) [2ng 



RPA 



(0)]'/',(7) 



where the superscript indicates that the RPA screening 
approximation has been used. In the limit that Ug 3> 
Tiimp, Eq. [71 reduces to Eq. [6l while in the opposite limit 
Ug — 7> 0, Eq. [3 reduces to results shown in Fig.[21 

Using the values for nimp and d determined from the 
bilayer data discussed above, and using Eq.[71for n^nis and 
Eq. [2l to calculate ^, we find theoretically (without any 



adjustable parameter) that 
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and 



^ = 8 nm, which should be compared to the experimental 
values extracted from the data reported in Ref. [l^. The 
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FIG. 6: Comparison of theoretical results with experimental 
data. Top panel shows the normalized correlation function 
A{r) = C(r)/C(0) for bilayer graphene. The circles are from 
the experimental data and solid curve is the theory for bi- 
layer graphene with d — 1 nm. The theory curve is insen- 
sitive to impurity concentration and doping away from the 
Dirac point. The error bars indicate single standard devi- 
ation uncertainties.— The small oscillation in the data over 
the monotonic decrease is a result of the finite size of the 
experimental image. Bottom panel is the normalized puddle 
correlation function in monolayer graphene at the Dirac point. 
Note the change in a;-axis scale from bilayer graphene in top 
panel. The solid curve is obtained from the self-consistent 
screening theory. The black squares are the results of a nu- 
merical mesoscopic density functional theory calculation for 
the ground-state properties of monolayer graphene^^ while 
the circles are experimental data taken from Deshpande et 
al.(Rei. Transport measurements on that same device 

set nimp = lO^^cm"^ which is the value used for the theory 
curves. The theory also uses d = 1 nm, which is the typical 
distance of the impurities from the graphene sheet extracted 
from transport measurements of graphene on Si02.— 



our theoretical calculations. This result should be con- 
trasted with bilayer graphene shown in Fig. [B^a) where 
the experimentally determined ^ — (3.68 ±0.03) nm, and 
the d = 1 nm RPA theory gives ^ = 3.4 nm. 

To further confirm our results, we compare our 
calculations to measurements made using scanning 
Coulomb blockade spectroscopy on a sample of mono- 
layer graphenCiii^ The circles in Fig. IHl^b) are experi- 
mental data for the normalized correlation and the solid 
line is the self-consistent theory discussed above us- 
ing riimp ~ lO^^cm^^ and d = 1 nm for monolayer 
graphene at the Dirac point. The impurity concentration 
and impurity distance were determined from transport 
measurementsiSiii^ and as such, no adjustable parameters 
were used in the calculation. 

In Fig. ini we also show (black squares) the results ex- 
tracted from a numerical mesoscopic density functional 
theory^ using the same parameters. The agreement be- 
tween the two calculations provides a posteriori justifi- 
cation for our assumption of a global screening function 
characterized by the density Ucb- 



V. CONCLUSIONS 

We conclude with the observation that our results re- 
quire only that the source of the disorder potential be 
uncorrelated charged impurities, and as such should ap- 
ply to graphene on other substrates. For example, re- 
cently graphene devices with hexagonal BN gate insula- 
tors have been fabricated showing transport properties 
similar to suspended graphene^'^ and larger puddles than 
on Si02 substrates .^^i^^ These observations are consis- 
tent with both a much smaller charged impurity density 



^imp 



on the BN substrate and with a larger distance d of 
the impurities from the graphene layer. Both these sce- 
narios are possible because the BN substrate is typically 
placed on top of the usual Si02 wafer which would have 
similar charged disorder to the samples we study here. 
We argue that an analysis similar to what we have per- 
formed here would be able to uniquely determine both d 
and riimp- Moreover, if a similar experiment is done with 
bilayer graphene on BN substrates, we predict that the 
puddle characteristics will not change much from what 
we find here with bilayer graphene on Si02. 
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